Demographics and reproductive biology of H. schistosus from bycatch on the west coast of India

Analysis for Manuscript

Introduction

This the code for analysis and supplementary figures for the manuscript. The manuscript is currently under review.

#Required libraries

require(tidyverse)
require(lubridate)
require(viridis)

# Importing data

snakes = read.csv("./Data/Sea-snakes_fish-dep_mastersheet_250420.csv")
embryos = read.csv("./Data/Sea snake_embryo data_mastersheet.csv")

snakes <- snakes%>%
  mutate(Year = as.factor(Year))%>%
  filter(Species == "Hydrophis curtus" | Species == "Hydrophis schistosus")%>%
  droplevels()

HS <- snakes%>%
  filter(Species == "Hydrophis schistosus")%>% # filtering out other snakes from bycatch data
  mutate(Year = as.factor(Year))

# function to calculate mcfadden's r2

mcfadden <- function(x){
  
  r2 <- 1 - (x$deviance/x$null.deviance)
  
}

# sampling years
years <- data.frame(Year = factor(c(2016:2019)))

Sample size

The number of snakes (N) sampled from bycatch of trawlers, shore seines and gillnetss.

snakes%>%
  group_by(Species)%>%
  count(name = "N")
Species N
Hydrophis curtus 282
Hydrophis schistosus 967

The number of snakes sampled in each year:

n.yr <- snakes%>%
  group_by(Species, Year)%>%
  count(name = "N")

n.yr%>%
  spread(Year, N)
Species 2016 2017 2018 2019
Hydrophis curtus 75 38 123 46
Hydrophis schistosus 178 215 521 53

The number of trips (n) and fishing effort (haul-hours) sampled for sea snakes in the current study.

snakes%>%
  filter(Gear.Type != "")%>%
  select(Date, Gear.Type, No..of.Hauls, Average.Haul.Duration..Hours., Tow.duration.hours)%>%
  distinct()%>%
  # Calculating tow duration
  mutate(Tow.duration.hours = case_when(!is.na(Tow.duration.hours) ~ Tow.duration.hours,
                                        is.na(Tow.duration.hours) ~ No..of.Hauls*Average.Haul.Duration..Hours.))%>%
  group_by(Gear.Type)%>%
  summarise(n = n(), # counting number of trips sampled
            haul.hours = sum(Tow.duration.hours, na.rm = T)) # total fishing effort
Gear.Type n haul.hours
GillNet 340 535.9633
Rampan 46 190.4500
Trawler 76 285.1600

Demographics of H. schistosus across time

Age structure

# calculating the mean SVL across years

age.yr <- snakes%>%
  group_by(Species, Year)%>%
  summarise(N = n(),
            mean = mean(Snout.to.Vent..cm., na.rm = T))

# SVL cut off for age classes

maturtity <- snakes%>%
  group_by(Species)%>%
  count()%>%
  mutate(juv = 35,
         adult = ifelse(Species == "Hydrophis curtus", 54, 65))

Distribution fo SVL across years

# plotting distribution of SVL across years

snakes%>%
  filter(Snout.to.Vent..cm. > 20)%>%
  ggplot(aes(Year, Snout.to.Vent..cm.))+
  geom_violin(fill = "grey")+
  geom_boxplot(width = 0.1)+
  geom_hline(data = maturtity, aes(yintercept = adult), linetype = "dashed")+
  geom_hline(data = maturtity, aes(yintercept = juv), linetype = "dotted")+
  #stat_summary(fun.data = "mean_sdl", geom = "pointrange", size = 1)+
  geom_label(data = age.yr, aes(Year, 10, label = N))+
  facet_wrap(~Species, ncol = 1, scales = "free_y")+
  labs(y = "Snout to vent length (cm)")+
  theme(strip.text = element_text(face = "italic"))

ggsave(last_plot(), filename = "./Figures/figure1.tiff", height = 6, width = 8)

We did not find any H. curtus neonates in our sampling. Majority of H. curtus were juveniles. H. schistosus populations consisted of mostly adults.

Linear model for change in SVL of each species across years:

library(car)

snakes%>%
  group_by(Species)%>%
  select(Year, Snout.to.Vent..cm.)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
          sumr = map(mod, broom::tidy),
         stat = map(mod, glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(c(Species:r.squared))
Species term estimate std.error statistic p.value r.squared
Hydrophis schistosus (Intercept) 85.598303 1.850695 46.2519793 0.0000000 0.0689121
Hydrophis schistosus Year2017 -13.996389 2.122893 -6.5930743 0.0000000 0.0689121
Hydrophis schistosus Year2018 -10.158177 1.974836 -5.1438083 0.0000003 0.0689121
Hydrophis schistosus Year2019 -1.732303 2.818894 -0.6145327 0.5390387 0.0689121
Hydrophis curtus (Intercept) 65.043404 1.511818 43.0233133 0.0000000 0.1470604
Hydrophis curtus Year2017 -5.212154 2.405287 -2.1669571 0.0312836 0.1470604
Hydrophis curtus Year2018 -8.833123 1.825449 -4.8388759 0.0000024 0.1470604
Hydrophis curtus Year2019 -12.998166 2.225335 -5.8409920 0.0000000 0.1470604

ANOVA to test for effect of years:

snakes%>%
  group_by(Species)%>%
  select(Year, Snout.to.Vent..cm.)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
          stat = map(mod, car::Anova))%>%
  select(stat)%>%
  unnest()
Species Sum Sq Df F value Pr(>F)
Hydrophis schistosus 13334.508 3 19.66265 0e+00
Hydrophis schistosus 180165.620 797 NA NA
Hydrophis curtus 4363.959 3 12.98867 1e-07
Hydrophis curtus 25310.652 226 NA NA

Plotting fitted values:

snakes%>%
  group_by(Species)%>%
  select(Year, Snout.to.Vent..cm.)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
         fit = map(mod, ~predict(., years, se.fit = T)),
         fit = map(fit, ~as.data.frame(.)))%>%
  select(fit)%>%
  unnest()%>%
  mutate(Year = factor(c(2016:2019)))%>%
  ggplot(aes(Year, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "SVL (cm)")+
  facet_wrap(~Species)+
  theme(strip.text = element_text(face = "italic"))

There was a significant change in SVL distribution of species across years.

Proportion of developmental classes across years

pop.yr <- snakes%>%
  filter(!is.na(Class))%>% # fix 2016 svl
  group_by(Species, Year)%>%
  # Counting number of individuals in each developmental class
  count(Class)%>%
  complete(Class, fill = list(n = 0))%>%
  # Total number sampled in each year
  mutate(N = sum(n))%>%
  group_by(Year, Species, Class)%>%
  # Proportion of each developmental class
  mutate(prop.age = n/N)

# Summary stats

pop.yr%>%
  group_by(Species, Class)%>%
  skimr::skim(prop.age)%>%
  skimr::yank("numeric")%>%
  select(c(Species, Class, mean, sd))

Variable type: numeric

Species Class mean sd
Hydrophis curtus Adult 0.43 0.31
Hydrophis curtus Juvenile 0.57 0.31
Hydrophis curtus Neonate 0.02 NA
Hydrophis schistosus Adult 0.88 0.09
Hydrophis schistosus Juvenile 0.14 0.07
Hydrophis schistosus Neonate 0.03 0.02

Testing change in proportion of developmental classes across years

Fitting a multinomial logistic model to determine change in age classes across years in each species.

require(nnet)
require(broom)

dev.yr <- snakes%>%
  select(Species, Year, Class)%>%
  filter(!is.na(Class))%>% # removing empty class 
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~multinom(Class ~ Year, data = .)),
         null = map(data, ~multinom(Class ~ 1, data = .)),
         smry = map(mod, ~broom::tidy(., exonentiate = T)),
         # Calculating mc fadden's r2
         r2 = map2_dbl(mod, null, ~ 1 - .x$value/.y$value))
dev.yr%>%
  select(smry, r2)%>%
  unnest()
Species y.level term estimate std.error statistic p.value r2
Hydrophis schistosus Juvenile (Intercept) -2.4450606 0.7380571 -3.3128339 0.0009236 0.0455695
Hydrophis schistosus Juvenile Year2017 1.1234010 0.7573074 1.4834147 0.1379643 0.0455695
Hydrophis schistosus Juvenile Year2018 0.5064641 0.7510767 0.6743174 0.5001096 0.0455695
Hydrophis schistosus Juvenile Year2019 -13.3315259 377.0158759 -0.0353606 0.9717922 0.0455695
Hydrophis schistosus Neonate (Intercept) -3.1352560 1.0212895 -3.0698995 0.0021413 0.0455695
Hydrophis schistosus Neonate Year2017 -14.5022162 526.2494814 -0.0275577 0.9780149 0.0455695
Hydrophis schistosus Neonate Year2018 -0.9347455 1.0900891 -0.8574946 0.3911716 0.0455695
Hydrophis schistosus Neonate Year2019 -8.7344420 53.4649902 -0.1633675 0.8702291 0.0455695
Hydrophis curtus Juvenile (Intercept) -1.9386218 1.0661369 -1.8183610 0.0690090 0.0724193
Hydrophis curtus Juvenile Year2017 2.7268690 1.1322946 2.4082681 0.0160284 0.0724193
Hydrophis curtus Juvenile Year2018 2.5058458 1.0853066 2.3088829 0.0209501 0.0724193
Hydrophis curtus Juvenile Year2019 3.5480223 1.1437099 3.1022049 0.0019208 0.0724193
Hydrophis curtus Neonate (Intercept) -11.1507510 99.7723289 -0.1117620 0.9110121 0.0724193
Hydrophis curtus Neonate Year2017 -2.2749071 278.6764080 -0.0081633 0.9934867 0.0724193
Hydrophis curtus Neonate Year2018 8.2062502 99.7749668 0.0822476 0.9344498 0.0724193
Hydrophis curtus Neonate Year2019 -3.4778373 576.2501607 -0.0060353 0.9951846 0.0724193

Plotting fitted estimates of proportions:

dev.yr%>%
  mutate(fit = map(mod, ~data.frame(predict(., years, "probs"))))%>%
  select(fit)%>%
  unnest()%>%
  mutate(Year = factor(c(2016:2019)))%>% # adding years variable
  left_join(n.yr, by = c("Species", "Year"))%>% # getting sample size in each year
  gather(c(Adult, Juvenile, Neonate), key = "Class", value = "p.hat")%>%
  mutate(se = sqrt(p.hat*(1-p.hat)/N))%>% # caluclating standard error
  # Plotting
  ggplot(aes(Year, p.hat, color = Species))+
  geom_pointrange(aes(ymin = p.hat-se, ymax = p.hat+se))+
  labs(y = "Estimated proportion")+
  facet_wrap(~Class)+
  theme(legend.text = element_text(face = "italic"))

Age structure of H. schistosus does not change significantly over a four year period from 2016 to 2018.

The proportion of H. curtus juveniles increases significantly between 2016 and 2018

Change in SVL distribution over months

month.svl <- HS%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(Month)%>%
  group_by(Month)%>%
  # Calculating mean SVL in each month
  summarise(mean.SVL = mean(Snout.to.Vent..cm., na.rm = T))

# Month of observed births

births <- data.frame(Species = "Hydrophis schistosus", Month = "April")
# plotting distribution of SVL across months

HS%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(nesting(Species), Month)%>%
  droplevels()%>%
  ggplot(aes(Month, Snout.to.Vent..cm.))+
  geom_violin(fill = "light grey")+
  #geom_point(data = month.svl, aes(x = Month, y = mean.SVL), size = 3)+
  geom_boxplot(width = 0.1)+
  geom_segment(data = births, aes(x = Month, xend = Month, y = 0 , yend = 10), #marking births
               arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
  geom_text(data = births, aes(x = Month, y = 20, label = paste("Observed \n births")))+
  geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
  geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
  geom_text(aes(x = "July", y = 80, label = "Monsoon ban"))+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  labs(y = "Snout to vent length (cm)")

Proportion of neonates increases in March to May. Birth and neonates were observed around the same time.

Proportion of neonates in each month of 2018

snakes%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  group_by(Species, Month)%>%
  summarise(prop.neonate = sum(Snout.to.Vent..cm. < 40)/n())%>%
  spread(Month, prop.neonate)
Species January February March April May October November December
Hydrophis curtus 0 0 0.50 NA NA 0 0.0769231 0
Hydrophis schistosus 0 0 0.02 0.05 0.0555556 0 0.0000000 0

Sex ratios

Proportion of females encountered by species

snakes%>%
  group_by(Species)%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  summarise(females = sum(Sex == "Female"),
            N = n(),
            P = females/N)
Species females N P
Hydrophis curtus 102 212 0.4811321
Hydrophis schistosus 304 618 0.4919094

Testing equality in sex ratios

# is proportion of females different from 0.5?

snakes%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  group_by(Species)%>%
  summarise(females = sum(Sex == "Female"),
            N = n())%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$females, .$N, p = 0.5)), # Z - test for proportion
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(-c(method, alternative))
Species estimate statistic p.value parameter conf.low conf.high
Hydrophis curtus 0.4811321 0.2311321 0.6306857 1 0.4125066 0.5504525
Hydrophis schistosus 0.4919094 0.1310680 0.7173273 1 0.4518628 0.5320580

Sex ratio is not different from 0.5; p = 0.71.

Testing shifts in sex ratios over years

snakes%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  group_by( Species, Year)%>%
  summarise(N = n(),
            females = sum(Sex == "Female"),
            prop.female = females/N)%>%
  ggplot(aes(Year, prop.female, fill = Species))+
  geom_col(width = 0.5, col = "black", position = position_dodge())

Logistic model to test for change in sex ratios of species across years:

snakes%>%
  select(Species, Year, Sex)%>%
  filter(Sex == "Male" | Sex == "Female")%>% # removing empty sex
  mutate(Sex = factor(Sex, levels = c("Male", "Female")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Sex ~ Year, family = binomial, data = .)),
         sumry = map(mod, ~tidy(.)),
         r2 = map_dbl(mod, ~mcfadden(.)))%>%
  select(sumry, r2)%>%
  unnest(sumry)
Species term estimate std.error statistic p.value r2
Hydrophis schistosus (Intercept) 0.5007753 0.2833779 1.7671643 0.0772007 0.0055997
Hydrophis schistosus Year2017 -0.6120009 0.4379020 -1.3975752 0.1622407 0.0055997
Hydrophis schistosus Year2018 -0.6052544 0.2977861 -2.0325137 0.0421017 0.0055997
Hydrophis schistosus Year2019 -0.3404326 0.4010216 -0.8489135 0.3959294 0.0055997
Hydrophis curtus (Intercept) 1.0185696 0.3235748 3.1478639 0.0016447 0.0593767
Hydrophis curtus Year2017 -1.3752445 0.5895404 -2.3327403 0.0196618 0.0593767
Hydrophis curtus Year2018 -1.4885732 0.3812161 -3.9048020 0.0000943 0.0593767
Hydrophis curtus Year2019 -1.2096248 0.4481189 -2.6993392 0.0069477 0.0593767

Plotting predicted values:

snakes%>%
  select(Species, Year, Sex)%>%
  filter(Sex == "Male" | Sex == "Female")%>% # removing empty sex
  mutate(Sex = factor(Sex, levels = c("Male", "Female")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Sex ~ Year, family = "binomial", data = .)),
         fit = map(mod, ~predict(., years, se.fit = T, type = "response")),
         fit = map(fit, ~data.frame(.)))%>%
  select(fit)%>%
  unnest(fit)%>%
  mutate(Year = factor(c(2016:2019)))%>%
  ggplot(aes(Year, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "P(Female|Year)")+
  facet_wrap(~Species)+
  theme(strip.text = element_text(face = "italic"))

The proportion of female H. curtus changed significantly across years. H. schistosus sex ratios remained constant across sampled years.

Mortality in bycatch

Overall mortality rate

snakes%>%
  group_by(Species)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Species n N prop.dead
Hydrophis curtus 120 282 0.4255319
Hydrophis schistosus 174 967 0.1799380

Testing difference in mortality across species

Logistic model:

snakes%>%
  filter(Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  nest()%>%
  mutate(test = map(data, ~glm(Condition ~ Species, family = "binomial", data = .)),
         sumr = map(test, ~broom::tidy(.)),
         r2 = map(test, ~mcfadden(.)))%>%
  select(sumr, r2)%>%
  unnest()
term estimate std.error statistic p.value r2
(Intercept) -0.2814125 0.1209241 -2.327182 0.0199556 0.0502489
SpeciesHydrophis schistosus -1.2315652 0.1470904 -8.372845 0.0000000 0.0502489

Fitted probability of mortality:

spp = data.frame(Species = c("Hydrophis schistosus", "Hydrophis curtus"))

snakes%>%
  filter(Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Condition ~ Species, family = "binomial", data = .)),
         fit = map(mod, ~predict(., spp, se.fit = T, type = "response")),
         fit = map(fit, ~data.frame(.)))%>%
  select(fit)%>%
  unnest(fit)%>%
  mutate(Species = factor(c("Hydrophis schsitosus", "Hydrophis curtus")))%>%
  ggplot(aes(Species, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "P(Dead|Species)")+
  theme(axis.text.x =  element_text(face = "italic"))

H. curtus had as significatly higher mortality rate in bycatch than H. schistosus

Mortality rate by age class

Are juveniles more vulnerable to bycatch mortality than adults?

snakes%>%
  filter(!is.na(Class))%>%
  group_by(Species, Class)%>%
  count(Condition.at.encounter..D.A.)%>%
  #calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Species Class n N prop.dead
Hydrophis curtus Adult 40 62 0.6451613
Hydrophis curtus Juvenile 41 125 0.3280000
Hydrophis schistosus Adult 138 648 0.2129630
Hydrophis schistosus Juvenile 5 105 0.0476190
Hydrophis schistosus Neonate 3 8 0.3750000

Testing difference in mortality rates across age classes within species

Logistic model:

snakes%>%
  filter(!is.na(Class),
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Class, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~glm(Condition ~ Class, family = "binomial", data = .)),
         sumr = map(test, ~broom::tidy(.)),
         r2 = map(test, ~mcfadden(.)))%>%
  select(sumr, r2)%>%
  unnest()
Species term estimate std.error statistic p.value r2
Hydrophis schistosus (Intercept) -1.3012573 0.0960144 -13.5527320 0.0000000 0.0300658
Hydrophis schistosus ClassJuvenile -1.6944749 0.4681739 -3.6193277 0.0002954 0.0300658
Hydrophis schistosus ClassNeonate 0.7904317 0.7365814 1.0731085 0.2832225 0.0300658
Hydrophis curtus (Intercept) 0.6931472 0.2738613 2.5310156 0.0113733 0.0823342
Hydrophis curtus ClassJuvenile -1.3984157 0.3338240 -4.1890813 0.0000280 0.0823342
Hydrophis curtus ClassNeonate -16.2592154 1029.1215009 -0.0157991 0.9873946 0.0823342

Fitted probability of mortality across age classes:

dev = data.frame(Class = c("Adult", "Juvenile", "Neonate"))

snakes%>%
  filter(!is.na(Class),
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Class, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Condition ~ Class, family = "binomial", data = .)),
         fit = map(mod, ~predict(., dev, se.fit = T, type = "response")),
         fit = map(fit, ~data.frame(.)))%>%
  select(fit)%>%
  unnest(fit)%>%
  mutate(Class = factor(c("Adult", "Juvenile", "Neonate")))%>%
  ggplot(aes(Class, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "P(Dead)")+
  facet_wrap(~Species)+
  theme(strip.text = element_text(face = "italic"))

H. curtus adults had the highest mortality rate overall. H. schistosus juveniles had the lowest.

Mortality rate by sex

snakes%>%
  filter(Sex != "")%>%
  group_by(Species, Sex)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Species Sex n N prop.dead
Hydrophis curtus Female 36 102 0.3529412
Hydrophis curtus Male 52 110 0.4727273
Hydrophis schistosus Female 62 304 0.2039474
Hydrophis schistosus Male 72 314 0.2292994

Testing difference in mortality rates across sexes withing species

Logistic model:

snakes%>%
  filter(Sex %in% c("Male", "Female"),
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Sex, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~glm(Condition ~ Sex, family = "binomial", data = .)),
         sumr = map(test, ~broom::tidy(.)),
         r2 = map(test, ~mcfadden(.)))%>%
  select(sumr, r2)%>%
  unnest()
Species term estimate std.error statistic p.value r2
Hydrophis schistosus (Intercept) -1.3535045 0.1424629 -9.5007493 0.0000000 0.0008565
Hydrophis schistosus SexMale 0.1453737 0.1957905 0.7424963 0.4577867 0.0008565
Hydrophis curtus (Intercept) -0.6061358 0.2071938 -2.9254529 0.0034396 0.0133815
Hydrophis curtus SexMale 0.5500463 0.2834464 1.9405655 0.0523110 0.0133815

Fitted probability of mortality across sexes:

sex = data.frame(Sex = c("Male", "Female"))

snakes%>%
  filter(Sex %in% c("Male", "Female"),
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Sex, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Condition ~ Sex, family = "binomial", data = .)),
         fit = map(mod, ~predict(., sex, se.fit = T, type = "response")),
         fit = map(fit, ~data.frame(.)))%>%
  select(fit)%>%
  unnest(fit)%>%
  mutate(Sex = c("Male", "Female"))%>%
  ggplot(aes(Sex, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "P(Dead)")+
  facet_wrap(~Species)+
  theme(strip.text = element_text(face = "italic"))

There was no difference in mortality rates across sexes in either species.

Mortality and female reproductive status

Are gravid females more susceptible to bycatch mortality than the rest of our sample?

HS%>%
  filter(Sex == "Female",
         Class == "Adult")%>%
  group_by(Gravid)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Gravid n N prop.dead
41 151 0.2715232
Yes 16 85 0.1882353

Testing thing difference in mortality rate of H. schistosus females by reproductive state

Logistic model:

HS%>%
  filter(Sex == "Female",
         Class == "Adult",
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Gravid, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")),
         Gravid = ifelse(Gravid == "", "No", Gravid))%>%
  nest()%>%
  mutate(test = map(data, ~glm(Condition ~ Gravid, family = "binomial", data = .)),
         sumr = map(test, ~broom::tidy(.)),
         r2 = map(test, ~mcfadden(.)))%>%
  select(sumr, r2)%>%
  unnest()
term estimate std.error statistic p.value r2
(Intercept) -0.9685592 0.1834379 -5.280039 0.0000001 0.0087809
GravidYes -0.4929586 0.3326292 -1.482006 0.1383386 0.0087809

Fitted probability of mortality by reproductive state:

gravid = data.frame(Gravid = c("Yes", "No"))

HS%>%
  filter(Sex == "Female",
         Class == "Adult",
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Gravid, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")),
         Gravid = ifelse(Gravid == "", "No", Gravid))%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Condition ~ Gravid, family = "binomial", data = .)),
         fit = map(mod, ~predict(., gravid, se.fit = T, type = "response")),
         fit = map(fit, ~data.frame(.)))%>%
  select(fit)%>%
  unnest(fit)%>%
  mutate(State = c("Gravid", "Not Gravid"))%>%
  ggplot(aes(State, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "P(Dead)")+
  theme(strip.text = element_text(face = "italic"))

Observed reproductive cycle of H. schistosus

Proportion of gravid females in the sample

# percentage of gravid females in sample

HS%>%
  count(Gravid)%>%
  mutate(N = sum(n))%>%
  mutate(prop.gravid = n/N)%>%
  filter(Gravid == "Yes")%>%
  select(-c(Gravid))
n N prop.gravid
88 967 0.0910031

Number of gravid females encountered per year

# checking the number of gravid females per year

HS%>%
  group_by(Year)%>%
  filter(Gravid == "Yes")%>%
  count(Gravid)%>%
  spread(Gravid, n)%>%
  rename(`n(Gravid Females)` = Yes)
Year n(Gravid Females)
2016 1
2018 75
2019 12

Monthly variation in proportion of gravid females

Proper sampling was only conducted in 2018/19 and hence only this period is used for analysis of reproductive cycles.

# calculating the proportion of gravid females per month

gravid <- HS%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  filter(Year == "2018" | Year == "2019")%>% # only for 2018/19
  group_by(Month)%>%
  summarise(N = n(),
            gravid = sum(Gravid == "Yes"),
            prop.gravid = gravid/N)

gravid
Month N gravid prop.gravid
January 102 12 0.1176471
February 134 34 0.2537313
March 129 24 0.1860465
April 54 7 0.1296296
May 36 1 0.0277778
October 16 0 0.0000000
November 69 3 0.0434783
December 34 6 0.1764706

Describing change in the relative proportions of gravid females across months and years of sampling.

# plotting prop gravid per month

gravid%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(Month, fill = list(prop.gravid = 0))%>%
  ggplot(aes(Month, prop.gravid))+
  geom_point(size = 3)+
  geom_path(aes(group = 1), size = 1, linetype = "dashed")+
  geom_segment(aes(x = "April", xend = "April", y = 0 , yend = 0.02), #marking births
               arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
  geom_text(aes(x = "April", y = 0.04, label = paste("Observed \n births")))+
  geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
  geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
  geom_text(aes(x = "July", y = 0.20, label = "Monsoon ban"))+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  labs(y = "Proportion of gravid females")

Pregnancy from Novemeber to May. Peak in Feb. Observed two live births in April.

Development of embryos and eggs

Sample size

embryos%>%
  summarise(n.mothers = length(unique(Field.Code)),
            n.embryos = length(unique(Embryo.Code)))
n.mothers n.embryos
29 235

Summary statistics of egg measurements

embryos%>%
  select(Egg.Length..mm.., Egg.Width..mm.., Egg.Weigth..g.., Snout.to.Vent..cm., Embryo.Weight..g.)%>%
  skimr::skim()%>%
  skimr::yank("numeric")%>%
  select(c(skim_variable, mean, sd))

Variable type: numeric

skim_variable mean sd
Egg.Length..mm.. 39.71 11.72
Egg.Width..mm.. 32.23 9.01
Egg.Weigth..g.. 14.30 5.45
Snout.to.Vent..cm. 17.72 5.24
Embryo.Weight..g. 4.83 6.28

Sex ratios within clutches

embryos%>%
  group_by(Field.Code)%>%
  filter(Sex != "")%>%
  summarise(prop.female = sum(Sex == "Female")/n())%>%
  skimr::skim(prop.female)%>%
  skimr::yank("numeric")%>%
  select(c(skim_variable, mean, sd))

Variable type: numeric

skim_variable mean sd
prop.female 0.53 0.31

Number of feamle embryos

embryos%>%
  filter(Sex != "")%>%
  summarise(females = sum(Sex == "Female"),
            N = n())
females N
86 166

Testing clutch sex ratio

broom::tidy(prop.test(86, 166, p = 0.5))%>% # Z - test
  select(estimate:conf.high)
estimate statistic p.value parameter conf.low conf.high
0.5180723 0.1506024 0.6979603 1 0.4395567 0.5957383

The sex ratio in clutches is not significantly different from 0.5.

Infertility rate

Percentage of eggs with no embryos present:

embryos%>%
  count(Embryo)%>%
  mutate(N = sum(n), rate = n*100/N)%>%
  select(-N)
Embryo n rate
Absent 5 2.12766
Present 230 97.87234

Matrotrophy

Do female H. schistosus expend energy in the development of embryos? or Are the eggs formed and yolk only provides nourishment?

Change in yolk weight with embryo weight

embryos%>%
  ggplot(aes(Embryo.Weight..g., Yolk.Weight..g.))+
  geom_point(size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  labs(x = "Embryo weight (g)", y = "Yolk weight (g)")

Linear model to test change in yolk weight with embryo development:

embryos%>%
  select(Yolk.Weight..g., Embryo.Weight..g.)%>%
  nest()%>%
  # Yolk weight is normally distributed
  mutate(mod = map(data, ~lm(Yolk.Weight..g. ~ Embryo.Weight..g., data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 13.0905565 0.2965229 44.14687 0 0.4071222
Embryo.Weight..g. -0.4526283 0.0450509 -10.04705 0 0.4071222

Yolk weight decreases by 0.45 g for every 1g increase in embryo weight.

Change in egg mass with embryo development

embryos%>%
  ggplot(aes(Embryo.Weight..g., Egg.Weigth..g..))+
  geom_point(aes(col = Yolk.Weight..g.), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_color_viridis(name = "Yolk Weight (g)")+
  labs(x = "Embryo Weight (g)", y = "Total Egg Weight (g)")

Linear model to test relation ship between total egg mass and embryo development:

embryos%>%
  select(Egg.Weigth..g.., Embryo.Weight..g.)%>%
  nest()%>%
  # Total egg mass is normally distributed
  mutate(mod = map(data, ~lm(Egg.Weigth..g.. ~ Embryo.Weight..g., data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 13.5441434 0.4114909 32.91481 0 0.2188866
Embryo.Weight..g. 0.3676922 0.0557912 6.59050 0 0.2188866

Total egg mass increases by 0.36 g for every 1g increase in embryo weight.

Note: Embryo weight is log-normally distributed.

Along with decrease of fat bodies of gravid females (Figure A1) as embryos develop indicates matrotrophic nutrition.

Reproductive effort

Does the amount of enery expended by female H. schistosus reduce with age?

# Cleaning reproductive effort data

re <- embryos%>%
  left_join(snakes, by = c("Date", "Field.Code", "Species"))%>%
  select(Date, Field.Code, Embryo.Code, Species, Egg.Length..mm..:Tail.Length..cm..x,
         Egg.Weigth..g..:Sex.x, Snout.to.Vent..cm..y:Tail.Length..cm..y, Weight..g.,
         -Body.Length..cm..x, - Body.Length..cm..y)%>%
  filter(Species == "Hydrophis schistosus")%>%
  rename(Embryo.SVL = Snout.to.Vent..cm..x,
         Embryo.TL = Tail.Length..cm..x,
         Embryo.Sex = Sex.x,
         Female.SVL = Snout.to.Vent..cm..y,
         Female.TL = Tail.Length..cm..y,
         Mother.Wt = Weight..g.
         )%>%
  # Removing erroneous data point
  filter(Female.SVL > 50)

Increase in clutch size with female age

re_clutch <- re%>%
  select(Field.Code, Egg.Weigth..g.., Mother.Wt, Female.SVL)%>%
  group_by(Field.Code)%>%
  summarise(clutch.size = n(),
            Clutch.wt = sum(Egg.Weigth..g..),
            Mother.wt = last(Mother.Wt),
            Female.SVL = last(Female.SVL))%>%
  mutate(rcm = Clutch.wt/(Mother.wt - Clutch.wt))

re_clutch%>%
  ggplot(aes(Female.SVL, clutch.size))+
  geom_point(size =3)+
  geom_smooth(method = 'glm', method.args = list(family = 'poisson'))+
  labs(x = "Female SVL (cm)", y = "Clutch Size")

Generalised linear model to test change in clutch size with female SVL:

Clutch size is poisson distributed with mean 8.1052632 and variance 10.3216374

re_clutch%>%
  select(Female.SVL, clutch.size)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(clutch.size ~ Female.SVL, data = ., family = "poisson")),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance),
         r.squared =  map_dbl(stat, ~mcfadden(.)))%>%
  select(sumr, stat,r.squared)%>%
  unnest()%>%
  select(c(term:p.value, r.squared))
term estimate std.error statistic p.value r.squared
(Intercept) -0.3764929 0.7174214 -0.5247862 0.5997318 0.560094
Female.SVL 0.0256926 0.0073002 3.5194451 0.0004325 0.560094

The number of eggs borne by females increase by 1.0253151 for every 1 cm increase in female SVL.

Change in overall reproductive effort with age

Beta regression to determine relationship between relative clutch mass and female SVL:

library(betareg)

re_clutch%>%
  select(Female.SVL, rcm)%>%
  nest()%>%
  mutate(mod = map(data, ~betareg(rcm ~ Female.SVL, data = .)),
         sumr = map(mod, ~tidy(., conf.int = T)),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(component:pseudo.r.squared)
component term estimate std.error statistic p.value conf.low conf.high pseudo.r.squared
mean (Intercept) 0.6463332 1.6515322 0.3913537 0.6955358 -2.5906105 3.8832769 0.0497975
mean Female.SVL -0.0144713 0.0172664 -0.8381195 0.4019636 -0.0483129 0.0193702 0.0497975
precision (phi) 12.8451280 4.7017084 2.7320129 0.0062949 3.6299489 22.0603071 0.0497975
re_mod <- betareg(rcm ~ Female.SVL, data = re_clutch)
re_clutch%>%
  mutate(fitted = predict(re_mod, re_clutch),
         q.low = predict(re_mod, re_clutch, "quantile", at = c(0.05)),# lower bound
         q.high = predict(re_mod, re_clutch, "quantile", at = c(0.95)))%>% # upper bound
  ggplot(aes(Female.SVL, rcm))+
  geom_point(aes(col = clutch.size), size = 3)+
  geom_smooth(aes(y = fitted, ymin = q.low, ymax = q.high), stat = "identity", linetype = "dashed")+
  scale_y_continuous(name = "Relative clutch mass")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Clutch size")+
  theme(legend.text = element_text(size = 10))

The overall reproductive effort does not change with female age.

Change in reproductive effort per embryo with female age

Does the effort per embryo change with female age?

# calculating reproductive effort per embruo

re_embryo <- re%>%
  select(Field.Code, Embryo.Code, Egg.Weigth..g..,  Mother.Wt, Female.SVL, Embryo.Sex, Embryo.SVL)%>%
  group_by(Field.Code)%>%
  mutate(clutch.wt = sum(Egg.Weigth..g..), 
         Female.wt = Mother.Wt - clutch.wt)%>%
  group_by(Field.Code, Embryo.Code)%>%
  summarise(inv = Egg.Weigth..g../Female.wt,# investment per embryo
            Female.SVL = last(Female.SVL),
            Embryo.Sex = Embryo.Sex,
            Embryo.SVL = Embryo.SVL)%>%
  ungroup()

Linear model to detemine relationship between female SVL and investment per embryo:

We have controlled for the effect of embryo development.

re_embryo%>%
  select(Female.SVL, inv, Embryo.SVL)%>%
  nest()%>%
  mutate(mod = map(data, ~betareg(inv ~ Female.SVL + Embryo.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr,stat)%>%
  unnest()%>%
  select(component:p.value, pseudo.r.squared)
component term estimate std.error statistic p.value pseudo.r.squared
mean (Intercept) -2.2798667 0.4894683 -4.657843 0.0000032 0.5627972
mean Female.SVL -0.0161182 0.0050085 -3.218169 0.0012901 0.5627972
mean Embryo.SVL 0.0491682 0.0077622 6.334311 0.0000000 0.5627972
precision (phi) 233.3499725 43.5079611 5.363386 0.0000001 0.5627972

Plotting relatonship between residuals and female SVL:

# using residuals to control for embryo development

emsvlinv <- betareg(inv ~ Embryo.SVL, data = re_embryo)

# Plotting

re_embryo%>%
  modelr::add_residuals(emsvlinv)%>%
  ggplot(aes(Female.SVL, resid))+
  geom_point(aes(col = Embryo.SVL), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_y_continuous(name = "Relative egg mass (residuals)")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Embryo SVL (cm)")

The relative egg mass (controlled for embryo development) reduces with female SVL.

Difference in reproductive effort for male and female embryos

Does female investment differ by sex of the embryo?

re_embryo%>%
  filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
  modelr::add_residuals(emsvlinv)%>%
  ggplot(aes(Female.SVL, resid, shape = Embryo.Sex))+
  geom_point(aes(col = Embryo.SVL), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_y_continuous(name = "Relative egg mass (residuals)")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Embryo SVL (cm)")+
  scale_shape_discrete(name = "Embryo Sex")

Linear model to test for the difference in female investment by embryo sex:

We controlled for the effect of embryo development.

re_embryo%>%
  filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
  select(Embryo.Sex, inv, Embryo.SVL)%>%
  nest()%>%
  mutate(mod = map(data, ~betareg(inv ~ Embryo.Sex + Embryo.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(component:p.value, pseudo.r.squared)
component term estimate std.error statistic p.value pseudo.r.squared
mean (Intercept) -3.7036341 0.1763377 -21.003076 0.0000000 0.4781387
mean Embryo.SexMale -0.1376065 0.0839749 -1.638663 0.1012835 0.4781387
mean Embryo.SVL 0.0481588 0.0083489 5.768291 0.0000000 0.4781387
precision (phi) 205.9146533 38.4135707 5.360466 0.0000001 0.4781387

Female investment does not differ significantly with embryo sex.